# -*- coding: utf-8 -*-
"""
PBAR_ProcessStandardWidthImagesVer2.py

4/7/2014, Added new boxes, added calculaion of value from average spectrum

Created on Thu Oct 10 13:59:03 2013

@author: jkwong
"""
import PBAR_Zspec, PBAR_Cargo
reload(PBAR_Zspec)
import numpy as np
import os,cPickle,copy
import matplotlib.pyplot as plt
#from scipy import interpolate
from matplotlib import cm
from scipy import ndimage


# Set useful plot variables
plotColors = ['r', 'b', 'g', 'm', 'c', 'y', 'k'] * 10
lineStyles = ['-', '-.', ':', '_', '|'] *  10
markerTypes = ['.', 'o', 'v', '^', '<', '>', '1', '2', '3', '4', 's', 'p', '*', 'h', 'H', '+', 'x', 'D', 'd']
markerTypes = markerTypes * 2
colorList = np.array(['r', 'b', 'g', 'm', 'c', 'y'])

(goodZspecMask, badZspecMask, goodZspecNameList, badZspecNameList, \
        goodZspecIndices, badZspecIndices) = PBAR_Zspec.ZspecDetectorLists2ndSet()

#dataPath = r'E:\PBAR\data4\BasicScansStandardWidth'
dataPath = r'E:\PBAR\data4\BasicScansStandardWidthPickle'

data0Path = r'E:\PBAR\data4\BasicScanCargo'

plotDir = r'E:\PBAR\data4\BasicScansPlots'
plot0Dir = r'E:\PBAR\data4\BasicScansPlots'

basepathSet2 = r'C:\Users\jkwong\Documents\Work\PBAR\data4\Mar-files'
plotSaveDir = r'E:\PBAR\data4\BasicScansPlots'

datasetDescription = PBAR_Zspec.ReadCargoDataDescriptionFile(r'E:\PBAR\data4\CargoSet2.txt')

#        cargoCountRange = np.array([0.01e8, 0.4e8]) # only consider parts of the zspec image that have rad below this value
cargoCountRange = np.array([0, 0.4e8]) # only consider parts of the zspec image that have rad below this value

acqTime = 1/60.
figureSize = (16, 10)

MINBIN = 1

savePlots = 1
closePlots = 1

# Load the coefficients and the polynomial of the discriminant band

## set 1        
#  wBest = np.array([-2.2806503 ,  2.26895887,  6.7898522 ])  # old 10/14/2013
#wBest = np.array([-2.07779891,  1.88095226,  4.65429156])
#pfitmean = np.array([ -3.45931698e-06,   1.80445624e-03,  -3.46939502e-01,  9.03548991e-01])
#pfitsigma = np.array([ -2.90915823e-07,   2.71344856e-04,  -6.81517865e-02, 6.54989380e+00])


# set 2 - load from file
setNum = 2
fullFilename = os.path.join(basepathSet2, 'zspec%dSet.dat' %setNum)
with open(fullFilename ,'rb') as fid:
    print('Reading %s' %fullFilename)
    temp = cPickle.load(fid)

polyMeanOrderExtendIndex = 22
polySigmaOrderExtendIndex = 13

pfitmean = temp['pfit']['mean']['PbALL']['extend'][polyMeanOrderExtendIndex]
pfitsigma = temp['pfit']['sigma']['PbALL']['extend'][polySigmaOrderExtendIndex]


#################################
## DEFINE FILE NAMES

energy = np.array([7.50E+00,1.05E+01,1.35E+01,1.65E+01,1.95E+01,2.25E+01,2.55E+01, \
    2.85E+01,3.15E+01,3.45E+01,3.75E+01,4.05E+01,4.35E+01,4.65E+01,4.95E+01,\
    5.25E+01,5.55E+01,5.85E+01,6.15E+01,6.45E+01,6.75E+01,7.05E+01,7.35E+01, \
    7.65E+01,7.95E+01,8.25E+01,8.55E+01,8.85E+01,9.15E+01,9.45E+01,9.75E+01,1.01E+02])

for datasetIndex, datasetDescript in enumerate(datasetDescription['scanID']):
    print(datasetIndex, datasetDescription['scanID'][datasetIndex])
#    if datasetIndex == 1:
#        break
#
#    if datasetIndex < 61:
#        continue
#
#    if datasetIndex > 64:
#        continue

    if datasetIndex == 63:
        continue


    acquisitionTime = acqTime
    
    ######################
    ## READ IN DATA
    
    # Read Zspec
#    filenameZspec = '%s-FDFC-All.npy' %datasetDescription['scanID'][datasetIndex]
#    fullFilenameZspec = os.path.join(data0Path, filenameZspec)
#    print('Loading %s' %fullFilenameZspec)
#    datZspec = np.load(fullFilenameZspec)
    
    filenameZspec = '%s-FDFC-All_SW.dat' %datasetDescription['scanID'][datasetIndex]
    fullFilenameZspec = os.path.join(dataPath, filenameZspec)
    print('Loading %s' %fullFilenameZspec)
    energy, datStandardWidth = PBAR_Zspec.ReadZspecBasicScanPickle(fullFilenameZspec)    
    
    # Read Cargo
#    filenameCargo = 'PBAR-%s.cargoimage' %datasetDescription['dataFile'][datasetIndex]
#    fullFilenameCargo = os.path.join(dataPath, filenameCargo)
#    print('Loading %s' %fullFilenameCargo)
##    datCargo= np.load(fullFilenameCargo)
#    (datCargo,bpp,formatt,flag,low1,high1,low2,high2) = PBAR_Cargo.ReadCargoImage(fullFilenameCargo)

    filenameCargo = 'PBAR-%s.cargoimageSW.dat' %datasetDescription['dataFile'][datasetIndex]
    fullFilenameCargo = os.path.join(dataPath, filenameCargo)
    print('Loading %s' %fullFilenameCargo)
    with open(fullFilenameCargo, 'rb') as fid:
        datCargoStandardWidth = cPickle.load(fid)

    # Read in marker files
    filenameMarker = filenameCargo.replace('cargoimageSW.dat', 'cargomarkerSW')
    fullFilenameMarker = fullFilenameCargo.replace('cargoimageSW.dat', 'cargomarkerSW')
    # some don't have marker files
    if os.path.exists(fullFilenameMarker):
#        markers = PBAR_Cargo.ReadCargoMarker(fullFilenameMarker)
        markerStandardWidth = PBAR_Cargo.ReadCargoMarker(fullFilenameMarker)
    else:
        markerStandardWidth = []
    
    #############################
    ##  CALCULATE FEATURES
    
    discrim = []
    temp = {}
    
    # interp version
#    discrim = PBAR_Zspec.CalculateFeaturesZspecBasicInterp(datStandardWidth, energy, MINBIN)   
    discrim = PBAR_Zspec.CalculateFeaturesZspecBasic(datStandardWidth, energy, MINBIN)

    # Calculate the dist, linear combination
    # make features array (num features X time slices X number of detectors)
#        featuresMatrix = np.zeros((3, discrim['binMean'].shape[0], discrim['binMean'].shape[1]))
#        featuresMatrix[0,:,:] = discrim['binMean']
#        featuresMatrix[1,:,:] = discrim['binSTD']
#        featuresMatrix[2,:,:] = discrim['multibin_20_ratio']
#        
#        discrim['dist0']  = featuresMatrix[0,:,:] *wBest[0] + featuresMatrix[1,:,:] *wBest[1] + featuresMatrix[2,:,:] *wBest[2]
    
    # 5/6/2014, Moved this to PBAR_Zspec.CalculateFeaturesZspecBasicInterp
    # new discriminant 4/28/2014        
#    featuresMatrix = np.zeros((2, discrim['binMean'].shape[0], discrim['binMean'].shape[1]))
#    featuresMatrix[0,:,:] = discrim['binSkew']
#    featuresMatrix[1,:,:] = discrim['binMean']
#
#    xbar_ = np.array([  1.28294017,  23.0544693 ])
#    scalings_ = np.array([[ 13.10009414],[ -0.03063372]])
#    discrim['dist0'] = (featuresMatrix[0,:,:] - xbar_[0]) * scalings_[0] + (featuresMatrix[0,:,:] - xbar_[1]) * scalings_[1]
    
#            # clean up bad values
#            cutt = np.isnan(dist0)
#            dist0[cutt] = 0.0

    # MAKE MASK
        
    # Cut in radiography image of standard size 
    cargoZspecStandardWidthMask = ( datCargoStandardWidth >= cargoCountRange[0] ) & \
        ( datCargoStandardWidth <= cargoCountRange[1] )
    
    # Clean up radiograph mask
    temp = cargoZspecStandardWidthMask
#    # do some stuff
#    temp = ndimage.binary_opening(temp)
#    cargoZspecStandardWidthCleanedMask = copy.copy(temp)
    cargoZspecStandardWidthCleanedMask = cargoZspecStandardWidthMask

    # output a npy file for each distriminant
    for (keyIndex, k) in enumerate(discrim.keys()):
        f = fullFilenameZspec.replace('.dat', '_%s.dat' %k)
        np.save(f, discrim[k])
        print('Wrote: %s' % f)

    markerValues = []
    for i, mark in enumerate(markerStandardWidth):
        markTemp = copy.deepcopy(mark)
        
        # CENTER POINT
        temp = {}
        x_center = mark['x']
        y_center = mark['y']
        
        temp['x'] = x_center
        temp['y'] = y_center
        temp['xmean'] = x_center
        temp['ymean'] = y_center
        
        # cycle through all features in the discrim dictionary
        discrimTemp = {}            
        for (keyIndex, key) in enumerate(discrim.keys()):
            discrimTemp[key] = discrim[key][x_center, y_center]
        discrimTemp['cargoCount'] = datCargoStandardWidth[x_center, y_center]
        temp['discrim'] = discrimTemp
        discrimTemp['cargoCount'] = datCargoStandardWidth[x_center, y_center]           
        
        markTemp['center'] = temp
        
        
        #  CENTER BLOCK
        # +/- 2 pixels in x direction
        # +/- 1 pixels in y direction
        # does not care about transmission mask
        temp = {}
        x_center = mark['x']
        y_center = mark['y']
        temp['x'] = np.arange(-2, 3) + x_center
        temp['y'] = np.arange(-1, 2) + y_center
        temp['xmean'] = x_center
        temp['ymean'] = y_center
        
        # make the bounds of the box
        x_range = np.array((temp['x'][0], temp['x'][-1])) # left < right
        y_range = np.array((temp['y'][0], temp['y'][-1])) # top  < bottom
        
        xarray = np.arange(discrim['dist0'].shape[0])
        yarray = np.arange(discrim['dist0'].shape[1])
        
        xcut = (xarray >= x_range[0]) & (xarray <= x_range[1])
        ycut = (yarray >= y_range[0]) & (yarray <= y_range[1])
        
        discrimTemp = {}
        discrimMeanTemp = {}
#        discrimMeanSpectrumTemp - {}
        
        for (keyIndex, key) in enumerate(discrim.keys()):
            discrimTemp[key] = discrim[key][:,goodZspecMask & ycut][xcut,:].flatten()  # cut the good pixels and flatten
            temp2 = discrim[key][:,goodZspecMask & ycut][xcut,:].flatten()
            cut = ~np.isnan(temp2) & ~np.isinf(temp2)
            discrimMeanTemp[key] = temp2[cut].mean()

        # calculate for mean spectrum
        # first reshape the spectrum
        tempSpectrum = datStandardWidth[:,goodZspecMask & ycut,:][xcut,:,:]  # include only good pixels
        tempReshaped = tempSpectrum.reshape((tempSpectrum.shape[0]*tempSpectrum.shape[1], tempSpectrum.shape[2]))
        tempReshapedMean = tempReshaped.mean(0)
        # calculate for mean spectrum
        discrimMeanSpectrumTemp = PBAR_Zspec.CalculateFeaturesZspecSingleSpectrum(tempReshapedMean, energy, MINBIN)
        
#        discrimTemp['cargoCount0'] = datCargoStandardWidth[:,ycut][xcut,:].flatten()
        discrimTemp['cargoCount'] = datCargoStandardWidth[:,goodZspecMask & ycut][xcut,:].flatten()
        
        print(discrimTemp['cargoCount'].shape)
#        discrimMeanTemp['cargoCount0'] = datCargoStandardWidth[:,ycut][xcut,:].flatten().mean()
        discrimMeanTemp['cargoCount'] = datCargoStandardWidth[:,goodZspecMask & ycut][xcut,:].flatten().mean()

#        discrimMeanSpectrumTemp['cargoCount0'] = datCargoStandardWidth[:,ycut][xcut,:].flatten().mean()
        discrimMeanSpectrumTemp['cargoCount'] = datCargoStandardWidth[:,goodZspecMask & ycut][xcut,:].flatten().mean()
        
        temp['discrim'] = discrimTemp
        temp['discrimMean'] = discrimMeanTemp
        temp['discrimMeanSpectrum'] = discrimMeanSpectrumTemp
        
        markTemp['centerBlock'] = temp
        del temp2, temp, discrimTemp, discrimMeanTemp, discrimMeanSpectrumTemp
        
        
        #  CENTER BLOCK TWO
        # +/- 4 pixels in x direction
        # Does not include adjacent detectors
        # does not include transmission mask
        temp = {}
        xBuffer = 4
        x_center = mark['x']
        y_center = mark['y']
        temp['x'] = np.arange(-xBuffer, xBuffer+1) + x_center
        temp['y'] = y_center
        temp['xmean'] = x_center
        temp['ymean'] = y_center
        
        # make the bounds of the box
        x_range = np.array((temp['x'][0], temp['x'][-1])) # left < right
        y_range = np.array((temp['y'], temp['y'])) # top  < bottom        
        
        xarray = np.arange(discrim['dist0'].shape[0])
        yarray = np.arange(discrim['dist0'].shape[1])
        
        # cut in time slices
        xcut = (xarray >= x_range[0]) & (xarray <= x_range[1])
        ycut = (yarray >= y_range[0]) & (yarray <= y_range[1])
         
        discrimTemp = {}
        discrimMeanTemp = {}
        
        for (keyIndex, key) in enumerate(discrim.keys()):
            discrimTemp[key] = discrim[key][:,goodZspecMask & ycut][xcut,:].flatten()  # cut the good pixels and flatten
            temp2 = discrim[key][:,goodZspecMask & ycut][xcut,:].flatten()
            cut = ~np.isnan(temp2) & ~np.isinf(temp2)
            discrimMeanTemp[key] = temp2[cut].mean()
            
        # calculate for mean spectrum
        # first reshape the spectrum
        tempSpectrum = datStandardWidth[:,goodZspecMask & ycut,:][xcut,:,:]  # include only good pixels
        tempReshaped = tempSpectrum.reshape((tempSpectrum.shape[0]*tempSpectrum.shape[1], tempSpectrum.shape[2]))
        tempReshapedMean = tempReshaped.mean(0) #average spectrum
        # calculate for mean spectrum
        discrimMeanSpectrumTemp = PBAR_Zspec.CalculateFeaturesZspecSingleSpectrum(tempReshapedMean, energy, MINBIN)            
        print('binMean shape')
        print(discrimTemp['binMean'].shape)
        
#        discrimTemp['cargoCount0'] = datCargoStandardWidth[xcut,:].flatten()
        discrimTemp['cargoCount'] = datCargoStandardWidth[:,goodZspecMask & ycut][xcut,:].flatten()
        print(discrimTemp['cargoCount'].shape)
#        discrimMeanTemp['cargoCount0'] = datCargoStandardWidth[xcut,:].flatten().mean()
        discrimMeanTemp['cargoCount'] = datCargoStandardWidth[:,goodZspecMask & ycut][xcut,:].flatten().mean()
#        discrimMeanSpectrumTemp['cargoCount0'] = datCargoStandardWidth[:,ycut][xcut,:].flatten().mean()
        discrimMeanSpectrumTemp['cargoCount'] = datCargoStandardWidth[:,goodZspecMask & ycut][xcut,:].flatten().mean()
        
        temp['discrim'] = discrimTemp
        temp['discrimMean'] = discrimMeanTemp
        temp['discrimMeanSpectrum'] = discrimMeanSpectrumTemp
        
        markTemp['centerBlock2'] = temp
        del temp2, temp, discrimTemp, discrimMeanTemp, discrimMeanSpectrumTemp
        
        
        #  MARKER BOX
        temp = {}
        temp['x'] = np.arange(mark['rec_left'], mark['rec_right']+1)
        temp['y'] = np.arange(mark['rec_top'], mark['rec_bottom']+1)
        temp['xmean'] = temp['x'].mean()
        temp['ymean'] = temp['y'].mean()
        
        # make the bounds of the box
        x_range = np.array((temp['x'][0], temp['x'][-1])) # left < right
        y_range = np.array((temp['y'][0], temp['y'][-1])) # top  < bottom
        
        xarray = np.arange(discrim['dist0'].shape[0])
        yarray = np.arange(discrim['dist0'].shape[1])
        
        xcut = (xarray >= x_range[0]) & (xarray <= x_range[1])
        ycut = (yarray >= y_range[0]) & (yarray <= y_range[1])
        
        discrimTemp = {}
        discrimMeanTemp = {}
        
        for (keyIndex, key) in enumerate(discrim.keys()):
            discrimTemp[key] = discrim[key][:,goodZspecMask & ycut][xcut,:].flatten()
            temp2 = discrim[key][:,goodZspecMask & ycut][xcut,:].flatten()
            cut = ~np.isnan(temp2) & ~np.isinf(temp2)
            discrimMeanTemp[key] = temp2[cut].mean()

        # calculate for mean spectrum
        # first reshape the spectrum
        tempSpectrum = datStandardWidth[:,goodZspecMask & ycut,:][xcut,:,:]  # include only good pixels
        tempReshaped = tempSpectrum.reshape((tempSpectrum.shape[0]*tempSpectrum.shape[1], tempSpectrum.shape[2]))  #flatten so that I can get mean
        tempReshapedMean = tempReshaped.mean(0)
        # calculate for mean spectrum
        discrimMeanSpectrumTemp = PBAR_Zspec.CalculateFeaturesZspecSingleSpectrum(tempReshapedMean, energy, MINBIN)            
            
#        discrimTemp['cargoCount0'] = datCargoStandardWidth[:,ycut][xcut,:].flatten()
        discrimTemp['cargoCount'] = datCargoStandardWidth[:,goodZspecMask & ycut][xcut,:].flatten()
#        discrimMeanTemp['cargoCount0'] = datCargoStandardWidth[:,ycut][xcut,:].flatten().mean()
        discrimMeanTemp['cargoCount'] = datCargoStandardWidth[:,goodZspecMask & ycut][xcut,:].flatten().mean()
#        discrimMeanSpectrumTemp['cargoCount0'] = datCargoStandardWidth[:,ycut][xcut,:].flatten().mean()
        discrimMeanSpectrumTemp['cargoCount'] = datCargoStandardWidth[:,goodZspecMask & ycut][xcut,:].flatten().mean()        
        
        temp['discrim'] = discrimTemp
        temp['discrimMean'] = discrimMeanTemp
        temp['discrimMeanSpectrum'] = discrimMeanSpectrumTemp
        
        markTemp['block'] = temp
        del temp2, temp, discrimTemp, discrimMeanTemp, discrimMeanSpectrumTemp
        
        
        #  MARKER BOX and is in cargo transmission mask
        temp = {}
        temp['x'] = np.arange(mark['rec_left'], mark['rec_right']+1)
        temp['y'] = np.arange(mark['rec_top'], mark['rec_bottom']+1)
        temp['xmean'] = temp['x'].mean()
        temp['ymean'] = temp['y'].mean()
        
        # make the bounds of the box
        x_range = np.array((temp['x'][0], temp['x'][-1])) # left < right
        y_range = np.array((temp['y'][0], temp['y'][-1])) # top  < bottom
        
        xarray = np.arange(discrim['dist0'].shape[0])
        yarray = np.arange(discrim['dist0'].shape[1])
        
        xcut = (xarray >= x_range[0]) & (xarray <= x_range[1])
        ycut = (yarray >= y_range[0]) & (yarray <= y_range[1])
        
        discrimTemp = {}
        discrimMeanTemp = {}
        
        maskCut = cargoZspecStandardWidthCleanedMask[:,goodZspecMask & ycut][xcut,:]    
        
        for (keyIndex, key) in enumerate(discrim.keys()):
            discrimTemp[key] = discrim[key][:,goodZspecMask & ycut][xcut,:][maskCut].flatten()
            temp2 = discrim[key][:,goodZspecMask & ycut][xcut,:][maskCut].flatten()
            cut = ~np.isnan(temp2) & ~np.isinf(temp2)
            discrimMeanTemp[key] = temp2[cut].mean()

        # calculate for mean spectrum
        # first reshape the spectrum
        tempReshaped = datStandardWidth[:,goodZspecMask & ycut,:][xcut,:,:][maskCut,:]  # include only good pixels
#        tempReshaped = temp.reshape((temp.shape[0]*temp.shape[1], temp.shape[2]))  #flatten so that I can get mean
        tempReshapedMean = tempReshaped.mean(0)
        # calculate for mean spectrum
        discrimMeanSpectrumTemp = PBAR_Zspec.CalculateFeaturesZspecSingleSpectrum(tempReshapedMean, energy, MINBIN)

#        discrimTemp['cargoCount0'] = datCargoStandardWidth[:,goodZspecMask & ycut][xcut,:].flatten() # includes counts outside of cargo mask
        discrimTemp['cargoCount'] = datCargoStandardWidth[:,goodZspecMask & ycut][xcut,:][maskCut].flatten()
#        discrimMeanTemp['cargoCount0'] = datCargoStandardWidth[:,goodZspecMask & ycut][xcut,:].flatten().mean()
        discrimMeanTemp['cargoCount'] = datCargoStandardWidth[:,goodZspecMask & ycut][xcut,:][maskCut].flatten().mean()
#        discrimMeanSpectrumTemp['cargoCount0'] = datCargoStandardWidth[:,goodZspecMask & ycut][xcut,:].flatten().mean()
        discrimMeanSpectrumTemp['cargoCount'] = datCargoStandardWidth[:,goodZspecMask & ycut][xcut,:][maskCut].flatten().mean()

        temp['discrim'] = discrimTemp
        temp['discrimMean'] = discrimMeanTemp
        temp['discrimMeanSpectrum'] = discrimMeanSpectrumTemp
        
        markTemp['blockMask'] = temp
        del temp2, temp, discrimTemp, discrimMeanTemp, discrimMeanSpectrumTemp      
        
        
#        #  MARKER DONUT BOX, SURROUNDING BOX
#        temp = {}
#        x_center = mark['x']
#        y_center = mark['y']
#        xBuffer = 8
#        yBuffer = 2
#        xWidth = 20
#        yWidth = 4
#
#        # inside box
#        x_range1 = np.array([mark['rec_left'] - xBuffer, mark['rec_right'] + xBuffer])
#        y_range1 = np.array([mark['rec_top'] - yBuffer, mark['rec_bottom'] + yBuffer])
#        
#        # other box
#        x_range2 = np.array([mark['rec_left'] - xBuffer - xWidth, mark['rec_right'] + xBuffer + xWidth])
#        y_range2 = np.array([mark['rec_top'] - yBuffer - yWidth, mark['rec_bottom'] + yBuffer + yWidth])
#        
#        #        temp['x'] = np.arange(mark['rec_left'], mark['rec_right']+1)
#        #        temp['y'] = np.arange(mark['rec_top'], mark['rec_bottom']+1)
#        #        temp['xmean'] = temp['x'].mean()
#        #        temp['ymean'] = temp['y'].mean()
#        #        
#        #        # make the bounds of the box
#        #        x_range = np.array((temp['x'][0], temp['x'][-1])) # left < right
#        #        y_range = np.array((temp['y'][0], temp['y'][-1])) # top  < bottom
#        #
#        xarray = np.arange(discrim['binMean'].shape[0])
#        yarray = np.arange(discrim['binMean'].shape[1])
#        
#        # donut mask
#        maskmask = np.zeros_like(discrim['binMean']) == 1  # all false
#        # set outer box to true and then inner box to false (exclude)
#        
#        # not sure why this is not working
#        #xcut = (xarray >= x_range2[0]) & (xarray <= x_range2[1])
#        #ycut = (yarray >= y_range2[0]) & (yarray <= y_range2[1])
#        #maskmask[xcut,:][:,ycut] = True
#        for xIndex in xarray:
#            for yIndex in yarray:
#                if (xIndex >= x_range2[0]) & (xIndex <= x_range2[1]):
#                    if (yIndex >= y_range2[0]) & (yIndex <= y_range2[1]):
#                        maskmask[xIndex, yIndex] = True 
#        
#        for xIndex in xarray:
#            for yIndex in yarray:
#                if (xIndex >= x_range1[0]) & (xIndex <= x_range1[1]):
#                    if (yIndex >= y_range1[0]) & (yIndex <= y_range1[1]):
#                        maskmask[xIndex, yIndex] = False
#        
#        discrimTemp = {}
#        discrimMeanTemp = {}
#        
##        maskCut = cargoZspecStandardWidthCleanedMask[:,goodZspecMask]
#        
#        maskmask_goodZspecMask = maskmask[:,goodZspecMask]
#        
#        for (keyIndex, key) in enumerate(discrim.keys()):
#            discrimTemp[key] = discrim[key][:,goodZspecMask][maskmask_goodZspecMask].flatten()
#            temp2 = discrim[key][:,goodZspecMask][maskmask_goodZspecMask].flatten()
#            cut = ~np.isnan(temp2) & ~np.isinf(temp2)
#            discrimMeanTemp[key] = temp2[cut].mean()
#        
#        # calculate for mean spectrum
#        # first reshape the spectrum
#        tempReshaped = datStandardWidth[:,goodZspecMask,:][maskmask_goodZspecMask,:]  # include only good pixels, this step already flattens
##        tempReshaped = temp.reshape((temp.shape[0]*temp.shape[1], temp.shape[2]))  #flatten so that I can get mean
#        tempReshapedMean = tempReshaped.mean(0)
#        # calculate for mean spectrum
#        discrimMeanSpectrumTemp = PBAR_Zspec.CalculateFeaturesZspecSingleSpectrum(tempReshapedMean, energy, MINBIN)
#
##        discrimTemp['cargoCount0'] = datCargoStandardWidth[maskmask].flatten()
#        discrimTemp['cargoCount'] = datCargoStandardWidth[:,goodZspecMask][maskmask_goodZspecMask].flatten()
##        discrimMeanTemp['cargoCount0'] = datCargoStandardWidth[maskmask].flatten().mean()
#        discrimMeanTemp['cargoCount'] = datCargoStandardWidth[:,goodZspecMask][maskmask_goodZspecMask].flatten().mean()        
#        discrimMeanSpectrumTemp['cargoCount'] = datCargoStandardWidth[:,goodZspecMask][maskmask_goodZspecMask].flatten().mean()           
#        
#        temp['discrim'] = discrimTemp
#        temp['discrimMean'] = discrimMeanTemp
#        temp['discrimMeanSpectrum'] = discrimMeanSpectrumTemp
#        
#        markTemp['dounut'] = temp
#        del temp2, temp, discrimTemp, discrimMeanTemp, discrimMeanSpectrumTemp
#        
        
        markerValues.append(markTemp)

    # convert to json and dump
    fullFilenameMarkerValues = fullFilenameMarker.replace('cargomarkerSW', 'cargomarkervaluesSW.dat')
    with open(fullFilenameMarkerValues, 'wb') as fid:
        cPickle.dump(markerValues, fid, 2)
        print('Write %s' %fullFilenameMarkerValues)
#        jsonText = json.dumps(markerValues)
#        with open(fullFilenameMarkerValues, 'wb') as fid:
#            fid.write(jsonText)
#            print('Write %s' %fullFilenameMarkerValues)
    
    ###################
    ## PLOTS 
    index = 0
    
    ## Standard Size Images ##
    
    ## Cargo standard size with markers
    # for zspec plot the total counts
    logIntensity = 1
    plotMarkers = True
    
    plt.figure(figsize = figureSize)
    plt.grid()
    
    if logIntensity:
        intensity = np.log(datCargoStandardWidth.T)
    else:
        intensity = datCargoStandardWidth.T
        
    plt.imshow(intensity, interpolation = 'nearest', aspect='auto', cmap = cm.Greys_r)
    plt.colorbar()
    plt.title(filenameCargo + ', Standard Size')
    
    if plotMarkers:
        for i, mark in enumerate(markerStandardWidth):
            x = mark['rec_left'], mark['rec_right'], mark['rec_right'], mark['rec_left'], mark['rec_left']
            y = mark['rec_bottom'], mark['rec_bottom'], mark['rec_top'], mark['rec_top'], mark['rec_bottom']
            plt.plot(x, y, 'r')
            plt.plot(mark['x'], mark['y'], 'xr')
            if 'left' in mark:
                plt.plot(mark['left']['x'], mark['left']['y'], 'gx')
            if 'right' in mark:
                plt.plot(mark['right']['x'], mark['right']['y'], 'go')
            plt.text(mark['x'], mark['y'],  mark['target'], fontsize = 12, color = 'r')

    plt.xlim((0, 2200))

    # Save image to file
    plotFullFilename = os.path.join(plotSaveDir, '%s_StandardSize.png' %filenameCargo)
    if savePlots:
        plt.savefig(plotFullFilename)
    if closePlots:
        plt.close('all')
    
    # Zspec image, counts with markers        
    plt.figure(figsize = figureSize)
    plt.grid()
    featureName1 = 'count'
    t = plt.imshow(discrim[featureName1].T, interpolation = 'nearest', aspect='auto', cmap = cm.Greys_r)
    plt.colorbar()
    if plotMarkers:
        for i, mark in enumerate(markerStandardWidth):
            x = mark['rec_left'], mark['rec_right'], mark['rec_right'], mark['rec_left'], mark['rec_left']
            y = mark['rec_bottom'], mark['rec_bottom'], mark['rec_top'], mark['rec_top'], mark['rec_bottom']
            plt.plot(x, y, 'r')
            plt.plot(mark['x'], mark['y'], 'xr')
            if 'left' in mark:
                plt.plot(mark['left']['x'], mark['left']['y'], 'gx')
            if 'right' in mark:
                plt.plot(mark['right']['x'], mark['right']['y'], 'go')
            plt.text(mark['x'], mark['y'],  mark['target'], fontsize = 12, color = 'r')
    plt.title(filenameZspec + ', Standard Size')
    plt.xlim((0, 2200))
            
    # save image to file
    plotFullFilename = os.path.join(plotSaveDir, '%s_%s_StandardSize.png' %(filenameZspec, featureName1))
    if savePlots:
        plt.savefig(plotFullFilename)
    if closePlots:
        plt.close('all')
            
    # Zspec dist image
    plt.figure(figsize = figureSize)
    plt.grid()
    featureName1 = 'dist0'
    t = plt.imshow(discrim[featureName1].T, interpolation = 'nearest', aspect='auto', cmap = cm.Greys_r)
    if plotMarkers:
        for i, mark in enumerate(markerStandardWidth):
            x = mark['rec_left'], mark['rec_right'], mark['rec_right'], mark['rec_left'], mark['rec_left']
            y = mark['rec_bottom'], mark['rec_bottom'], mark['rec_top'], mark['rec_top'], mark['rec_bottom']
            plt.plot(x, y, 'r')
            plt.plot(mark['x'], mark['y'], 'xr')
            if 'left' in mark:
                plt.plot(mark['left']['x'], mark['left']['y'], 'gx')
            if 'right' in mark:
                plt.plot(mark['right']['x'], mark['right']['y'], 'go')
            plt.text(mark['x'], mark['y'],  mark['target'], fontsize = 12, color = 'r')
    plt.colorbar()
    plt.title('%s, %s ' %(filenameZspec, featureName1))
    plt.xlim((0, 2200))
    
    plotFullFilename = os.path.join(plotSaveDir, '%s_%s_StandardSize.png' %(filenameZspec, featureName1))
    if savePlots:
        plt.savefig(plotFullFilename)
    
    
    # Zspec dist image with mask
    plt.figure(figsize = figureSize)
    plt.grid()
    featureName1 = 'dist0'
    im = copy.copy(discrim[featureName1].T)
    minValue = im.min()
    im[~cargoZspecStandardWidthCleanedMask.T] = minValue
    t = plt.imshow(im, interpolation = 'nearest', aspect='auto', cmap = cm.Greys_r)
    if plotMarkers:
        for i, mark in enumerate(markerStandardWidth):
            x = mark['rec_left'], mark['rec_right'], mark['rec_right'], mark['rec_left'], mark['rec_left']
            y = mark['rec_bottom'], mark['rec_bottom'], mark['rec_top'], mark['rec_top'], mark['rec_bottom']
            plt.plot(x, y, 'r')
            plt.plot(mark['x'], mark['y'], 'xr')
            if 'left' in mark:
                plt.plot(mark['left']['x'], mark['left']['y'], 'gx')
            if 'right' in mark:
                plt.plot(mark['right']['x'], mark['right']['y'], 'go')
            plt.text(mark['x'], mark['y'],  mark['target'], fontsize = 12, color = 'r')
    
    plt.colorbar()
#    plt.title('Linear Combination')
    plt.title('%s, %s Linear Combination with Mask' %(filenameZspec, featureName1))
    plt.xlim((0, 2200))
    plt.show()

    plotFullFilename = os.path.join(plotSaveDir, '%s_%s_Mask_StandardSize.png' %(filenameZspec, featureName1))
    if savePlots:
        plt.savefig(plotFullFilename)
    if closePlots:
        plt.close('all')

    # Cargo mask
    plt.figure(figsize = figureSize)
    plt.grid()
    t = plt.imshow(cargoZspecStandardWidthCleanedMask.T, interpolation = 'nearest', aspect='auto', cmap = cm.Greys_r)
    if plotMarkers:
        for i, mark in enumerate(markerStandardWidth):
            x = mark['rec_left'], mark['rec_right'], mark['rec_right'], mark['rec_left'], mark['rec_left']
            y = mark['rec_bottom'], mark['rec_bottom'], mark['rec_top'], mark['rec_top'], mark['rec_bottom']
            plt.plot(x, y, 'r')
            plt.plot(mark['x'], mark['y'], 'xr')
            if 'left' in mark:
                plt.plot(mark['left']['x'], mark['left']['y'], 'gx')
            if 'right' in mark:
                plt.plot(mark['right']['x'], mark['right']['y'], 'go')
            plt.text(mark['x'], mark['y'],  mark['target'], fontsize = 12, color = 'r')
    
    plt.colorbar()
    plt.title('%s, Cargo Mask: %d, %d' %(filenameCargo, cargoCountRange[0], cargoCountRange[1]))
    plt.xlim((0, 2200))
    plt.show()

    plotFullFilename = os.path.join(plotSaveDir, '%s_%s_CargoMask_StandardSize.png' %(filenameZspec, featureName1))
    if savePlots:
        plt.savefig(plotFullFilename)
    if closePlots:
        plt.close('all')



    ### Scatter Plots ###

    #  All data points and mask, and those in the boxes
    plt.figure(figsize = figureSize)
    plt.grid()
    index = 0
    feature1 = 'count'
    feature2 = 'dist0'
    # x, y for the whole image
    x1 = discrim[feature1][:,goodZspecIndices].flatten()
    y1 = discrim[feature2][:,goodZspecIndices].flatten()
    
    # reduce the mask to good detectors
    cut = cargoZspecStandardWidthCleanedMask[:, goodZspecIndices]
    # x values in mask, good detectors        
    x2 = discrim[feature1][:,goodZspecIndices][cut].flatten()
    # y values in mask, good detectors        
    y2 = discrim[feature2][:,goodZspecIndices][cut].flatten()
    
    plt.plot(x1, y1, '.k', alpha = 0.15, markersize = 10, label = 'All Points')
    plt.plot(x2, y2, '.b', alpha = 0.15, markersize = 10, label = 'In Mask')
    
    # plot the points in the markers
    for i, mark in enumerate(markerStandardWidth):
        x_range = np.array((mark['rec_left'], mark['rec_right'])) # left < right
        y_range = np.array((mark['rec_top'], mark['rec_bottom'])) # top  < bottom
        xarray = np.arange(discrim[feature1].shape[0])
        yarray = np.arange(discrim[feature1].shape[1])
        xcut = (xarray > x_range[0]) & (xarray < x_range[1])
        ycut = (yarray > y_range[0]) & (yarray < y_range[1])
        
        x = discrim[feature1][:,goodZspecMask & ycut][xcut,:].flatten()
        y = discrim[feature2][:,goodZspecMask & ycut][xcut,:].flatten()
        
    #            plt.plot(x, y, '.', markersize = 10, alpha  = 0.3, marker = markerTypes[i/6], color = colorList[i%6], label = mark['target'])
        # high density stuff
        try:
            if mark['target'][0] == 'S':
                plt.plot(x, y, 'db', markersize = 8, alpha  = 0.3, label = mark['target'])
            elif mark['target'][0] == 'W':
                plt.plot(x, y, 'vb', markersize = 8, alpha  = 0.3, label = mark['target'])
            elif mark['target'][0] == 'D':
                plt.plot(x, y, 'sb', markersize = 8, alpha  = 0.3, label = mark['target'])
            elif mark['target'][0] == 'P':
                plt.plot(x, y, '*b', markersize = 8, alpha  = 0.3, label = mark['target'])
            elif mark['target'][0] == 'F':# low density stuff
                plt.plot(x, y, 'or', markersize = 8, alpha  = 0.3, label = mark['target'])
            else:
                plt.plot(x, y, 'om', markersize = 8, alpha  = 0.3, label = mark['target'])            
        except:
            plt.plot(x, y, 'om', markersize = 8, alpha  = 0.3, label = mark['target'])
        # high density stuff
    #    xx, yy = markerValues[i]['centerBlock']['discrimMean'][feature1], markerValues[i]['centerBlock']['discrimMean'][feature2]
        xx, yy = markerValues[i]['center']['discrim'][feature1], markerValues[i]['center']['discrim'][feature2]
    
        try:
            if markerValues[i]['target'][0] == 'S':
                plt.plot(xx, yy, 'db', markersize = 18, alpha  = 1.0, label = mark['target'] +' Center')
            elif mark['target'][0] == 'W':
                plt.plot(xx, yy, 'vb', markersize = 18, alpha  = 1.0, label = mark['target']+' Center')
            elif mark['target'][0] == 'D':
                plt.plot(xx, yy, 'sb', markersize = 18, alpha  = 1.0, label = mark['target']+' Center')
            elif mark['target'][0] == 'P':
                plt.plot(xx, yy, '*b', markersize = 18, alpha  = 1.0, label = mark['target']+' Center')
            elif mark['target'][0] == 'F':# low density stuff
                plt.plot(xx, yy, 'or', markersize = 18, alpha  = 1.0, label = mark['target']+' Center')
            else:
                plt.plot(xx, yy, 'om', markersize = 18, alpha  = 1.0, label = mark['target']+' Center')
        except:
            plt.plot(x, y, 'om', markersize = 8, alpha  = 0.3, label = mark['target'])
    
    plt.title(filenameZspec + ', Linear Combination Vs Counts')
    plt.legend(prop={'size':8})
    plt.xlabel(feature1)
    plt.ylabel(feature2)
    plt.xlim((0, 500))
    plt.ylim((-50, 50))
    plt.show()

    plotFullFilename = os.path.join(plotSaveDir, '%s_%s_%s_AllPoints_Mask_Markers.png' %(filenameZspec, featureName1, feature2))
    if savePlots:
        plt.savefig(plotFullFilename)
    if closePlots:
        plt.close('all')

    #  Data points in mask only 
    removeOffset = 0
    
    plt.figure(figsize = figureSize)
    plt.grid()
    feature1 = 'count'
    feature2 = 'dist0'
    # reduce the mask to good detectors
    cut = cargoZspecStandardWidthCleanedMask[:, goodZspecIndices]
    # x values in mask, good detectors        
    x2 = discrim[feature1][:,goodZspecIndices][cut].flatten()
    # y values in mask, good detectors        
    y2 = discrim[feature2][:,goodZspecIndices][cut].flatten()
    if removeOffset:
        y2 = y2 - np.polyval(pfitmean, x2)
    plt.plot(x2, y2, '.k', alpha = 0.15, markersize = 10, label = 'In Mask')
    
    # plot the points in the markers
    for i, mark in enumerate(markerStandardWidth):
        x_range = np.array((mark['rec_left'], mark['rec_right'])) # left < right
        y_range = np.array((mark['rec_top'], mark['rec_bottom'])) # top  < bottom
        xarray = np.arange(discrim[feature1].shape[0])
        yarray = np.arange(discrim[feature1].shape[1])
        xcut = (xarray > x_range[0]) & (xarray < x_range[1])
        ycut = (yarray > y_range[0]) & (yarray < y_range[1])
        x = discrim[feature1][:,goodZspecMask & ycut][xcut,:].flatten()
        y = discrim[feature2][:,goodZspecMask & ycut][xcut,:].flatten()
    #            plt.plot(x, y, '.', markersize = 15, alpha  = 0.3, marker = markerTypes[i/6], color = colorList[i%6], label = mark['target'])
        if removeOffset:
            y = y - np.polyval(pfitmean, x)
            
        try:
            # high density stuff
            if mark['target'][0] == 'S':
                plt.plot(x, y, 'db', markersize = 8, alpha  = 0.3, label = mark['target'])
            elif mark['target'][0] == 'W':
                plt.plot(x, y, 'vb', markersize = 8, alpha  = 0.3, label = mark['target'])
            elif mark['target'][0] == 'D':
                plt.plot(x, y, 'sb', markersize = 8, alpha  = 0.3, label = mark['target'])
            elif mark['target'][0] == 'P':
                plt.plot(x, y, '*b', markersize = 8, alpha  = 0.3, label = mark['target'])
            elif mark['target'][0] == 'F':# low density stuff
                plt.plot(x, y, 'or', markersize = 8, alpha  = 0.3, label = mark['target'])
            else:
                plt.plot(x, y, 'om', markersize = 8, alpha  = 0.3, label = mark['target'])
        except:
            plt.plot(x, y, 'om', markersize = 8, alpha  = 0.3, label = mark['target'])
        
        # high density stuff
        #    xx, yy = markerValues[i]['centerBlock']['discrimMean'][feature1], markerValues[i]['centerBlock']['discrimMean'][feature2]
        xx, yy = markerValues[i]['center']['discrim'][feature1], markerValues[i]['center']['discrim'][feature2]
        
        if removeOffset:
            yy = yy - np.polyval(pfitmean, xx)
        try:
            if markerValues[i]['target'][0] == 'S':
                plt.plot(xx, yy, 'db', markersize = 18, alpha  = 1.0, label = mark['target'] +' Center')
            elif mark['target'][0] == 'W':
                plt.plot(xx, yy, 'vb', markersize = 18, alpha  = 1.0, label = mark['target']+' Center')
            elif mark['target'][0] == 'D':
                plt.plot(xx, yy, 'sb', markersize = 18, alpha  = 1.0, label = mark['target']+' Center')
            elif mark['target'][0] == 'P':
                plt.plot(xx, yy, '*b', markersize = 18, alpha  = 1.0, label = mark['target']+' Center')
            elif mark['target'][0] == 'F':# low density stuff
                plt.plot(xx, yy, 'or', markersize = 18, alpha  = 1.0, label = mark['target']+' Center')
            else:
                plt.plot(xx, yy, 'om', markersize = 18, alpha  = 1.0, label = mark['target']+' Center')
        except:
            plt.plot(xx, yy, 'om', markersize = 18, alpha  = 1.0, label = mark['target']+' Center')
            
    plt.title(filenameZspec + ', Linear Combination Vs Counts, Special Points Not Necessarily in Mask')
    plt.legend(prop={'size':8})
    plt.xlabel(feature1)
    plt.ylabel(feature2)
    plt.xlim((0, 500))
    plt.ylim((-50, 50))
    plt.show()
    
    plotFullFilename = os.path.join(plotSaveDir, '%s_%s_%s_Mask_Markers.png' %(filenameZspec, featureName1, feature2))
    if savePlots:
        plt.savefig(plotFullFilename)
    if closePlots:
        plt.close('all')

    
    #  Data points in mask, box points can contain points outside mask
    removeOffset = 1
    
    plt.figure(figsize = figureSize)
    plt.grid()
    feature1 = 'count'        
    feature2 = 'dist0'
    # reduce the mask to good detectors
    cut = cargoZspecStandardWidthCleanedMask[:, goodZspecIndices]
    # x values in mask, good detectors        
    x2 = discrim[feature1][:,goodZspecIndices][cut].flatten()
    # y values in mask, good detectors        
    y2 = discrim[feature2][:,goodZspecIndices][cut].flatten()
    if removeOffset:
        y2 = y2 - np.polyval(pfitmean, x2)
    plt.plot(x2, y2, '.k', alpha = 0.15, markersize = 10, label = 'In Mask')
    
    # plot the points in the markers
    for i, mark in enumerate(markerStandardWidth):
        x_range = np.array((mark['rec_left'], mark['rec_right'])) # left < right
        y_range = np.array((mark['rec_top'], mark['rec_bottom'])) # top  < bottom
        xarray = np.arange(discrim[feature1].shape[0])
        yarray = np.arange(discrim[feature1].shape[1])
        xcut = (xarray > x_range[0]) & (xarray < x_range[1])
        ycut = (yarray > y_range[0]) & (yarray < y_range[1])
        x = discrim[feature1][:,goodZspecMask & ycut][xcut,:].flatten()
        y = discrim[feature2][:,goodZspecMask & ycut][xcut,:].flatten()
    #            plt.plot(x, y, '.', markersize = 15, alpha  = 0.3, marker = markerTypes[i/6], color = colorList[i%6], label = mark['target'])
        if removeOffset:
            y = y - np.polyval(pfitmean, x)
        try:
            # high density stuff
            if mark['target'][0] == 'S':
                plt.plot(x, y, 'db', markersize = 8, alpha  = 0.3, label = mark['target'])
            elif mark['target'][0] == 'W':
                plt.plot(x, y, 'vb', markersize = 8, alpha  = 0.3, label = mark['target'])
            elif mark['target'][0] == 'D':
                plt.plot(x, y, 'sb', markersize = 8, alpha  = 0.3, label = mark['target'])
            elif mark['target'][0] == 'P':
                plt.plot(x, y, '*b', markersize = 8, alpha  = 0.3, label = mark['target'])
            elif mark['target'][0] == 'F':# low density stuff
                plt.plot(x, y, 'or', markersize = 8, alpha  = 0.3, label = mark['target'])
            else:
                plt.plot(x, y, 'om', markersize = 8, alpha  = 0.3, label = mark['target'])
        except:
            plt.plot(x, y, 'om', markersize = 8, alpha  = 0.3, label = mark['target'])
        
        # high density stuff
        #    xx, yy = markerValues[i]['centerBlock']['discrimMean'][feature1], markerValues[i]['centerBlock']['discrimMean'][feature2]
        xx, yy = markerValues[i]['center']['discrim'][feature1], markerValues[i]['center']['discrim'][feature2]
        
        if removeOffset:
            yy = yy - np.polyval(pfitmean, xx)
        try:
            if markerValues[i]['target'][0] == 'S':
                plt.plot(xx, yy, 'db', markersize = 18, alpha  = 1.0, label = mark['target'] +' Center')
            elif mark['target'][0] == 'W':
                plt.plot(xx, yy, 'vb', markersize = 18, alpha  = 1.0, label = mark['target']+' Center')
            elif mark['target'][0] == 'D':
                plt.plot(xx, yy, 'sb', markersize = 18, alpha  = 1.0, label = mark['target']+' Center')
            elif mark['target'][0] == 'P':
                plt.plot(xx, yy, '*b', markersize = 18, alpha  = 1.0, label = mark['target']+' Center')
            elif mark['target'][0] == 'F':# low density stuff
                plt.plot(xx, yy, 'or', markersize = 18, alpha  = 1.0, label = mark['target']+' Center')
            else:
                plt.plot(xx, yy, 'om', markersize = 18, alpha  = 1.0, label = mark['target']+' Center')
        except:
            plt.plot(xx, yy, 'om', markersize = 18, alpha  = 1.0, label = mark['target']+' Center')
    plt.title(filenameZspec + ', Linear Combination Vs Counts, Special Points Not Necessarily in Mask')
    plt.legend(prop={'size':8})
    plt.xlabel(feature1)
    plt.ylabel(feature2)
    plt.xlim((0, 500))
    plt.ylim((-50, 50))
    plt.show()
    
    plotFullFilename = os.path.join(plotSaveDir, '%s_%s_%s_Mask_Markers_OffsetRemoved.png' %(filenameZspec, featureName1, feature2))
    if savePlots:
        plt.savefig(plotFullFilename)
    if closePlots:
        plt.close('all')
    
    
    #  Data points in mask, box points can contain points outside mask
    removeOffset = 1
    
    plt.figure(figsize = figureSize)
    plt.grid()
    feature1 = 'count'        
    feature2 = 'dist0'
    # reduce the mask to good detectors
    cut = cargoZspecStandardWidthCleanedMask[:, goodZspecIndices]
    # x values in mask, good detectors        
    x2 = discrim[feature1][:,goodZspecIndices][cut].flatten()
    # y values in mask, good detectors        
    y2 = discrim[feature2][:,goodZspecIndices][cut].flatten()
    if removeOffset:
        y2 = y2 - np.polyval(pfitmean, x2)
    plt.plot(x2, y2, '.k', alpha = 0.15, markersize = 10, label = 'In Mask')
    
    # plot the points in the markers
    for i, mark in enumerate(markerStandardWidth):
        x_range = np.array((mark['rec_left'], mark['rec_right'])) # left < right
        y_range = np.array((mark['rec_top'], mark['rec_bottom'])) # top  < bottom
        xarray = np.arange(discrim[feature1].shape[0])
        yarray = np.arange(discrim[feature1].shape[1])
        xcut = (xarray > x_range[0]) & (xarray < x_range[1])
        ycut = (yarray > y_range[0]) & (yarray < y_range[1])

        # only includ points in Box AND Mask
        cut = cargoZspecStandardWidthCleanedMask[:,goodZspecMask & ycut][xcut,:]

        x = discrim[feature1][:,goodZspecMask & ycut][xcut,:][cut].flatten()
        y = discrim[feature2][:,goodZspecMask & ycut][xcut,:][cut].flatten()
        
    #            plt.plot(x, y, '.', markersize = 15, alpha  = 0.3, marker = markerTypes[i/6], color = colorList[i%6], label = mark['target'])
        if removeOffset:
            y = y - np.polyval(pfitmean, x)
        try:
            # high density stuff
            if mark['target'][0] == 'S':
                plt.plot(x, y, 'db', markersize = 8, alpha  = 0.3, label = mark['target'])
            elif mark['target'][0] == 'W':
                plt.plot(x, y, 'vb', markersize = 8, alpha  = 0.3, label = mark['target'])
            elif mark['target'][0] == 'D':
                plt.plot(x, y, 'sb', markersize = 8, alpha  = 0.3, label = mark['target'])
            elif mark['target'][0] == 'P':
                plt.plot(x, y, '*b', markersize = 8, alpha  = 0.3, label = mark['target'])
            elif mark['target'][0] == 'F':# low density stuff
                plt.plot(x, y, 'or', markersize = 8, alpha  = 0.3, label = mark['target'])
            else:
                plt.plot(x, y, 'om', markersize = 8, alpha  = 0.3, label = mark['target'])
        except:
            plt.plot(x, y, 'om', markersize = 8, alpha  = 0.3, label = mark['target'])
        
        # high density stuff
        #    xx, yy = markerValues[i]['centerBlock']['discrimMean'][feature1], markerValues[i]['centerBlock']['discrimMean'][feature2]
        xx, yy = markerValues[i]['center']['discrim'][feature1], markerValues[i]['center']['discrim'][feature2]
        
        if removeOffset:
            yy = yy - np.polyval(pfitmean, xx)
        try:
            if markerValues[i]['target'][0] == 'S':
                plt.plot(xx, yy, 'db', markersize = 18, alpha  = 1.0, label = mark['target'] +' Center')
            elif mark['target'][0] == 'W':
                plt.plot(xx, yy, 'vb', markersize = 18, alpha  = 1.0, label = mark['target']+' Center')
            elif mark['target'][0] == 'D':
                plt.plot(xx, yy, 'sb', markersize = 18, alpha  = 1.0, label = mark['target']+' Center')
            elif mark['target'][0] == 'P':
                plt.plot(xx, yy, '*b', markersize = 18, alpha  = 1.0, label = mark['target']+' Center')
            elif mark['target'][0] == 'F':# low density stuff
                plt.plot(xx, yy, 'or', markersize = 18, alpha  = 1.0, label = mark['target']+' Center')
            else:
                plt.plot(xx, yy, 'om', markersize = 18, alpha  = 1.0, label = mark['target']+' Center')
        except:
            plt.plot(xx, yy, 'om', markersize = 18, alpha  = 1.0, label = mark['target']+' Center')
    plt.title(filenameZspec + ', Linear Combination Vs Counts, Special Points Not Necessarily in Mask')
    plt.legend(prop={'size':8})
    plt.xlabel(feature1)
    plt.ylabel(feature2)
    plt.xlim((0, 500))
    plt.ylim((-50, 50))
    plt.show()
    
    plotFullFilename = os.path.join(plotSaveDir, '%s_%s_%s_MaskOnly_Markers_OffsetRemoved.png' %(filenameZspec, featureName1, feature2))
    if savePlots:
        plt.savefig(plotFullFilename)  
    if closePlots:
        plt.close('all')
    